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Abstract. We analyse protostellar mass accretion rates M from numerical models of star formation based on gravoturbulent 
fragmentation, considering a large number of different environments. To within one order of magnitude, M a Mj /rg with Mj 
being the mean thermal Jeans mass and Tjf the corresponding free-fall time. However, mass accretion rates are highly time- 
variant, with a sharp peak shortly after the formation of the protostellar core. We present an empirical exponential fit formula to 
describe the time evolution of the mass accretion and discuss the resulting fit parameters. There is a positive correlation between 
the peak accretion rate and the final mass of the protostar. We also investigate the relation of M with the turbulent flow velocity 
as well as with the driving wavenumbers in different environments. We then compare our results with other theoretical models 
of star formation and with observational data. 
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1. Introduction 

Stars are born in dense cores of interstellar molecular clouds. 
Despite recent observational and theoretical progress, the ini- 
tial conditions and physical processes that determine the for- 
mation of stars are still not fully understood. 

In the so-called "standard theory of star formation" (Shu 
et al. 11987) stars are formed by the inside-out collapse of a 
singular isothermal sphere that is initially in quasistatic equi- 
librium, supported against gravity by magnetic and thermal 
pressure and evolves only due to slow ambipolar diffusion pro- 
cesses. This model predicts protostellar mass accretion rates 
that are constant with time and only depend on the isothermal 
sound speed (Shu '1977'). This hypothesis, however, has been 
challenged from several sides (see Larson 2003 or Mac Low & 
Klessen 120041 for a summary). It only is applicable to isolated, 
single stars, while it is known that the majority of stars form 
in small aggregates or large clusters (Adams & Myers 2001, 
Lada & Lada l200'3l . Furthermore, there is both observational 
evidence (Crutche/T^, Andre et al.'SUnU' Bourke et ak^ MHl 
and theoretical reasoning (e.g. Nakano 1 998) showing that 
most observed cloud cores do not have magnetic fields strong 
enough to support against gravitational collapse. Similarly, the 
long lifetimes implied by the quasi-static phase of evolution 
in the model are diflicult to reconcile, e.g., with observational 



statistics of cloud cores (Taylor et al. 1996 Lee & Myers|T999l 
Visser et al. .2002) and with chemical age considerations (van 
Dishoeck & Blake lT^ Langer et al. lITO . 



Molecular clouds appear to actually be transient objects 
with lifetimes of a few million years that form and dissolve in 
the larger scale turbulent flow of the Galactic disc (Ballesteros- 
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Hartmann 2003; Vazquez-Semadeni et al. l2004l l. Observations 
of self-similar structure in molecular clouds (e.g. Mac Low & 
Ossenkopf 2000 , Ossenkopf & Mac Low 2002 1 indicate that in- 
terstellar turbulence is driven on scales substantially larger than 
the clouds themselves. These large-scale turbulent flows com- 
press and cool gas. At sufficiently high densities atomic gas is 
then quickly converted into molecular form (Hollenbach et al. 
il971,) . These same flows will continue to drive the turbulent 
motions observed within the newly formed cloud. Some com- 
bination of turbulent flow, free expansion at the sound speed of 
the cloud, and dissociating radiation from internal star forma- 
tion will then be responsible for their destruction on a timescale 
of 5-10 Myr The most likely source of such large-scale inter- 
stellar turbulence in the Milky Way is the combined energy and 
momentum input from supernovae explosions. They appear to 
overwhelm all other possibilities. In the outer reaches of the 
Galaxy and in low surface brightness galaxies, on the other 
hand, the situation is not so clear, with magnetorotational or 
gravitational instabilities looking most likely to drive the ob- 
served flows (Mac Low .2002: Mac Low & Klessen .2004i) . 

Modern star formation theory, therefore, considers super- 
sonic interstellar turbulence as controlling agent for stellar 
birth, rather than mediation by magnetic fields (Mac Low 
& Klessen '2004). This turbulence typically carries suffi- 
cient energy to counterbalance gravity on global scales. On 
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small scales, however, it may actually provoke localised col- 
lapse (Hunter & Fleck IT9821 Elmegreen [T9931 Padoan 1995 
Ballesteros-Paredes et al. ll999al Klessen et al. '2000i Padoan & 
Nordlund .1999. .20021 1 . This apparent paradox can be resolved 
when considering that supersonic turbulence establishes a com- 
plex network of interacting shocks, where converging flows 
generate regions of enhanced density. The system becomes 
highly filamentary, with elongated structures being caused ei- 
ther by shear motions or by compression at the intersection of 
shocked layers of gas. At some locations the density enhance- 
ment can be sufficiently strong for gravitational instability to 
set in. The stability criteria for filaments and sheets have been 
derived and discussed in the context of star formation, e.g., by 
Larson STW5\ , Lubow & Pringle STW^ , and Clarke (|l999l. 
However, the same random flow that creates density enhance- 
ments may disperse them again. For local collapse to result in 
stellar birth, it must progress fast enough for the region to 'de- 
couple' from the flow. 

The efficiency of protostellar core formation, the growth 
rates and final masses of the protostars, and essentially all other 
properties of nascent star clusters then depend on the intricate 
interplay between gravity on the one hand side and the turbu- 
lent velocity field in the cloud on the other. The star forma- 
tion rate is regulated not just at the scale of individual star- 
forming cores through ambipolar diffusion balancing magne- 
tostatic support, but rather at all scales (Elmegreen l2002> . via 
the dynamical processes that determine whether regions of gas 
become unstable to prompt gravitational collapse. The pres- 
ence of magnetic fields does not alter that picture significantly 
(Mac Low et al.l9W, Stone et al. [T998l Padoan & Nordlund 
1 1999 Heitsch et al. 2001b i. In particular, it cannot prevent the 
decay of interstellar turbulence. 

Clusters of stars build up in molecular cloud regions where 
self-gravity overwhelms turbulence, either because such re- 
gions are compressed by a large-scale shock, or because in- 
terstellar turbulence is not replenished and decays on short 
timescales. Then, many gas clumps become gravitationally un- 
stable synchronously to go into collapse. If the number density 
is high, contracting protostellar cores interact and may merge 
to produce new cores which now contain multiple protostars. 
Close encounters drastically alter the trajectories of the proto- 
stars, thus changing their mass accretion rates. This has impor- 
tant consequences for the final stellar mass spectrum (Bonnell 
et al. \TW7\ Klessen & Burkert HUU^ EUUTl Bonnell et al. 
I2001a'"2001b Klessen 2001b"i Bate et al.'2002^. 

Inefficient, isolated star formation will occur in regions 
which are supported by turbulence carrying most of its energy 
on very small scales. This requires an unrealistically large num- 
ber of driving sources and appears at odds with the measured 
velocity structure in molecular clouds which in almost all cases 
is dominated by large-scale modes (Mac Low & Ossenkopf 
I2OOOI Ossenkopf & Mac Low 2002 1. 

In this paper we extend the analysis of protostellar mass ac- 
cretion rates from a single case (Klessen"2001a l to a large series 
of numerical models of turbulent molecular cloud fragmenta- 
tion, which essentially cover the entire spectrum of observed 
star-forming regions, ranging from inefficient and isolated star 
formation to the fast and efficient build-up of stellar clusters. 



These calculations, their numerical realisation, and the adopted 
parameters are described in Section|2] In Section|3lwe discuss 
our findings. We investigate the mass growth history of all pro- 
tostars in our set of models and present a simple analytic fit 
formula for the accretion rate M. We discuss our study in rela- 
tion with previous analyses and observational data in Sections 
|3and|5] respectively. Finally, in Section |6l we summarise our 
results. 

2. The models 

To adequately describe the fragmentation of turbulent, self- 
gravitating gas clouds, and the resulting formation and mass 
growth of protostars, it is prerequisite to resolve the dynami- 
cal evolution of collapsing cores over several orders of mag- 
nitude in density. Due to the stochastic nature of supersonic 
turbulence, it is not known in advance where and when this 
local collapse occurs. Hence, SPH (smoothed particle hydro- 
dynamics) is used to solve the equations of hydrodynamics. It 
is a Lagrangian method, where the fluid is represented by an 
ensemble of particles and flow quantities are obtained by av- 
eraging over an appropriate subset of the SPH particles (Benz 
^990 Monaghan ll992> . The method is able to resolve large 
density contrasts as particles are free to move and so naturally 
the particle concentration increases in high-density regions. 

We use the same smoothing procedure for gravity and pres- 
sure forces. This is one requirement to prevent artificial frag- 
mentation (Bate & Burkert .1997) . Because it is computation- 
ally prohibitive to treat the cloud as a whole, we concentrate on 
subregions within the cloud and adopt periodic boundary con- 
ditions (Klessen 1997 ) . Once the central region of a collapsing 
protostellar core exceeds a density contrast of ~ 10^^, it is re- 
placed by a "sink" particle (Bate et al. 1995 1, which has the 
ability to accrete gas from its surrounding while at the same 
time keeping track of mass and linear and angular momentum. 
By adequately replacing high-density cores with sink particles 
we can follow the dynamical evolution of the system over many 
free-fall times. 

The suite of models consists of two globally unstable mod- 
els that contract from Gaussian initial conditions without tur- 
bulence (for details see Klessen & Burkert 120001 120011 and of 
22 models where turbulence is maintained with constant rms 
Mach numbers Ai, in the range 0.1 < Al < 10. We distin- 
guish between turbulence that carries its energy mostly on large 
scales, at wavenumbers 1 < < 2, on intermediate scales, i.e. 
3 < k <4, and on small scales with 7 < < 8. The coiTespond- 
ing wavelengths are £ - L/k, where L is the total size of the 
computed volume. The models are labelled mnemonically as 
MAIkA;, with rms Mach number M and wavenumber k, while 
Gl and G2 denote the two Gaussian runs. The main parameters 
are summarised in Table ^ To have well defined environmen- 
tal conditions given by Ai and k, M is required to be constant 
throughout the evolution. However, turbulent energy dissipates 
rapidly, roughly on a free-fall timescale (Mac Low et al. 1 19981 
Stone et al. 1998. Padoan & Nordlund 1999 1. We therefore ap- 
ply a non-local driving scheme that inserts energy at a given 
rate and at a given scale k. We use Gaussian random fluctua- 
tions in velocity. This is appealing because Gaussian fields are 
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Table 1. Overview of our models (See text for details) 
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" number of particles in the simulation 
* SPH resolution limit 

^ fraction of the total mass that has been accreted by the end of the simulation 

total number of formed protostars 
'' number of protostars that can be fitted by Eq. jsj 
^ mean deviation of the fits, calculated from Eq. J6} 



fully determined by their power distribution in Fourier space. 
We define a cartesian mesh with 64^ cells, and for each three- 
dimensional wave number k we randomly select an amplitude 
from a Gaussian distribution around unity and a phase between 
zero and 27r. We then transform the resulting field back into real 
space to get a "kick-velocity" in each cell. Its amplitude is de- 
termined by solving a quadratic equation such to keep A1 con- 
stant (Mac Low [D^ Klessen et aL EUnOt . The "kick-velocity" 
is then simply added to the speed of each SPH particle located 
in the cell. We adopted this method for mathematical simplic- 
ity. In reality, the situation is far more complex. Still, our mod- 
els of large-scale driven clouds contain many features of molec- 
ular clouds in supernovae driven turbulence (e.g. Ballesteros- 
Paredes & Mac Low lMHl Mac Low et al. l2003> . Conversely, 
our models of small-scale turbulence bear certain resemblance 
to energy input on small scales provided by protostellar feed- 
back via outflows and winds. 

Our models neglect the influence of magnetic fields, be- 
cause their presence cannot halt the decay of turbulence 
(Mac Low et al. [T998l Stone et al. [T998l Padoan & Nordlund 
I1999j and does not significantly alter the efficiency of local 



collapse for driven turbulence (Heitsch et al. 12001 al l. More im- 
portantly, we do not self-consistently consider feedback effects 
from the star formation process itself (like bipolar outflows, 
stellar winds, or ionising radiation from new-born O or B stars). 
Our analysis of protostellar mass accretion rates solely focuses 
on the interplay between turbulence and self-gravity only. This 
is also the case in the Shu (.1977. 1 theory of isothermal collapse. 
Hence, our findings can be directly compared to the "standard 
theory of star formation". 



The models are computed in normalised units using an 
isothermal equation of state. Scaled to physical units we adopt 
a temperature of 1 1.3 K corresponding to a sound speed - 
0.2 km s"', and we use a mean density of n(H2) = 10^ cm"-*, 
which is typical for star-forming molecular cloud regions (e.g. 
in p Ophiuchi, see Motte et al. 199? I. The total mass con- 
tained in the computed volume in the two Gaussian models 
is 220 M0 and the size of the cube is 0.34 pc. This corre- 
sponds to 220 thermal Jeans masses. The turbulent models 
have a mass of 120 M© within a volume of (0.28 pc)^, equiv- 
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alent to 120 thermal Jeans masses'. The mean thermal Jeans 
mass in all models is thus (Mj) - 1 Mq, the global free-fall 
timescale is fff = 10^ yr, and the simulations cover a density 
range from n(H2) ~ 100 cm"^ in the lowest density regions 
to n{H2) ~ 10^ cm"^ where collapsing protostellar cores are 
identified and converted into "sink" particles in the code. This 
coincides in time with the formation of the central protostar 
to within ~ 10^ yr (Wuchterl & Klessen l2001> . The resolution 
limit for each model, requiring that the local Jeans mass is al- 
ways resolved by at least 100 gas particles (Bate & Burkert 
nWl\ . is given in Col. 5 of Tabled 

In the subsequent protostellar phase of evolution, we de- 
termine accretion rates M by measuring the amount of mass 
as function of time that falls into a control volume defined by 
each "sink" particle. Its diameter is fixed to 560 AU. Entering 
gas particles pass through several tests to check if they remain 
bound to the "sink" particle (Bate et al. I1995> before they are 
considered accreted. As all gas particles have the same mass 
and as accretion events occur at random times, the resulting 
accretion rates are mass-binned and we smooth over a few con- 
secutive accretion events to get a description of the time evolu- 
tion of M. We cannot resolve the evolution in the interior of the 
control volume. Because of angular momentum conservation 
most of the matter that falls in will assemble in a protostellar 
disc. There it is transported inwards by viscous and possibly 
gravitational torques. The latter will be provided by spiral den- 
sity waves that develop when the disc becomes too massive, 
which happens when mass is loaded onto the disc faster than it 
is removed by viscous transport alone. Altogether, the disc will 
not prevent or delay material from accreting onto the protostar 
for long. It acts as a buffer and smoothes eventual accretion 
spikes. For the mass range considered here also feedback ef- 
fects are too weak to halt or delay accretion. With typical disc 
sizes of order of several hundred AU, the control volume there- 
fore fully encloses both, star and disc, and the measured core 
accretion rates are good estimates for the actual stellar accre- 
tion rates. Deviations may be expected only if the protostellar 
core forms a binary star, where the infalling mass must then be 
distributed between two stars, or if very high-angular momen- 
tum material is accreted, where a certain mass fraction may end 
up in a circumbinary disc and not accrete onto a star at all. 

In the prestellar phase, i.e. before the central protostar 
forms, we determine the accretion history by computing the 
change of mass inside the control volume centered on the 
SPH particle that turns into a "sink" during the later evolution. 
Turbulent compression leads to mass growth, i.e. M > 0, while 
expansion will result in mass loss and M < 0. Appreciable 
mass growth, however, is only achieved when gravity takes 
over and the region goes into collapse. 



We use a spherical definition of the Jeans mass, Mj = 



1/2 



4/37Tp{A}/2y, with density p and Jeans length As = and 
where G and K are the gravitational and the gas constant. The mean 
Jeans mass (Afj) is then determined from the average density in the 
system (p). 



3. Discussion 

3.1. First approximation to M 

The entire process of molecular cloud collapse and build-up 
of the stellar cluster lasts several global free-fall times (% - 
10^ yr). Likewise, the accretion process of a protostellar core 
takes place on a timescale of a few fff, comparable to most 
other models of star formation. 

A simple approximation to the accretion rate can be 
achieved by dividing the local Jeans mass by the local dynam- 
ical timescale: 



M K Mi Its 
By substituting 



M> = —Po 



1/2 



3/2 ./2 



(1) 



(2) 



where po denotes the initial density, T the temperature and Cs 
the isothermal sound speed, and 



Tff = 



37T 

'32Gpo 



Eq. Q can be written as 



Ml 
M« — = 

Tff 



32 Tp- cj _ , 



— — ^ =5.4^, 
3 6 G G 



(3) 



(4) 



depending only on the isothermal sound speed (or tempera- 
ture). For a sound speed Cs = 0.2 kms"' we obtain M - 
10"^ Mq yr"' . This is higher than the accretion rate for the col- 
lapse of a singular isothermal sphere: M = 0.975 c^/G (Shu 
119771 . However, the accretion rates in our models vary with 
time. Typical peak accretion rates are roughly in the range 
(3 - 50) cl/G or 5 x 10"^ to lO^^^ Mq yr"' . The values exceed 
the approximated value Mj/xff due to external compression in 
the turbulent flow. 



3.2. Time-varying mass accretion rates 

We analyse the full mass growth history of all protostellar cores 
in our models and we find that mass accretion rates from gravo- 
turbulent fragmentation are highly time-variable. Several ex- 
amples of the accretion rate M are displayed in Fig. [2 plot- 
ted versus time (left panel) and the ratio of accreted to final 
mass (right panel), respectively. The maximum accretion rate 
is reached rather rapidly and is then followed by a somewhat 
slower decline. In some cases this decline is inteiTupted by one 
or more secondary peaks. As shown above, the maximum ac- 
cretion rate is significantly higher than the constant rate pre- 
dicted by the classical isothermal collapse model (plotted as 
dotted line in Fig.^, but it falls below that value in later stages. 
Due to the dynamical interaction and competition between pro- 
tostellar cores, the mass accretion rates of cores in a dense clus- 
ter are different from those of isolated cores. In the first stage a 
core accretes local gas from its immediate vicinity. Once the lo- 
cal reservoir is depleted, the core may accrete fresh gas stream- 
ing in from farther away or by encounters with non-collapsed 
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Fig. 1. Mass accretion rates of nine randomly selected protostellar cores of three different models. Left panel: M versus time for a 
Gaussian collapse (G2; upper row), a turbulent model driven on a large scale (M6k2a; middle row), and a turbulent model driven 
on a small scale (M6k8b; lower row). The thin line represents the actual simulation, the thick line the fit as described in the text. 
The deviation cr as given by Eq. (|6} is indicated for each object. The dotted line shows the constant accretion rate that would be 
expected from the classical Shu (I1977> scenario. The dashed line stands for the assumed transition from Class to Class I. The 
right panel shows the same protostellar cores as on the left side plotted versus the ratio of accreted to final mass. The final masses 
(in Mq) are also given. 



gas clumps (see discussion in Klessen & Burkert l2000> . This 
results in secondary accretion peaks that are also visible in the 
right panel of Fig. [2 where one would expect a single bump in 
the case of an isolated core. For example, the central graph of 
the right panel of Fig.^nicely shows that this particular pro- 
tostar accretes only about half of the final mass from its direct 
environment (first bump), while the rest stems from later accre- 
tion events. 

The transition phase between Class and Class I protostars 
is believed to take place when about half of the final mass has 
been accumulated (Andre et al. 2000 1. This time is indicated by 
the dashed line in Fig. ^ Typically it takes place during or at 
the end of the peak accretion phase. It determines the lifetime 
of Class objects, which will be discussed below. 

We define a mean accretion rate (M) by averaging M in the 
mass range 0.1 < M/Mend ^ 0.8, with Mend being the final 
mass of the protostar. This phase typically lasts only a few 10'* 
years. This is short compared to the full accretion history. The 
bulk of stellar material is therefore accumulated in the short 
time interval while the system is close to maximum accretion. 

In Fig. 12 we plot the mean accretion rates versus final star 
mass Mend and versus time of core formation ffoi-m. respectively, 
for the same models as in Fig. ^ Not surprisingly, (M) in- 
creases with increasing stellar mass, and decreases when the 
core forms later, although this second coiTelation is not that 
clear. In other words, more massive stars have higher mass ac- 
cretion rates and start to form first. They can grow large, be- 
cause on average they form in the high-density regions of the 



cluster centre where they are able to maintain relatively high 
accretion rates over a long time as more and more gas falls in 
from the cluster outskirts. 



3.3. An empirical fit formula for M 

One of our aims is to find a simple-to-use fit formula to ap- 
proximate the accretion process. The protostellar mass growth 
history in our models can be fitted empirically in the lin-log 
diagram by the function 



log M(f) = log M() - t e 

T 



-tiT 



(5) 



with time t and the fit parameters log Mo and t. This 
holds for the following conditions: We shift the ordinate by 
AlogM/(M0yr"') = +1 and we consider accretion when 
logM > -7. The zero point of the timescale is determined 
once the accretion rate exceeds log M - -1 . The fitted curves 
are plotted as thick lines in Fig. ^ Note that the ordinate dis- 
plays the original values without the applied shift. If there are 
secondary accretion peaks, they are typically ignored and lev- 
elled out by the routine. The overall quality of the fit can be 
estimated by the standard deviation 



1 " 



(6) 



where M(f) is the actual value of M at the time t from our sim- 
ulation, while Mfit(f) denotes M calculated using Eq. (|5} for 
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timescales on which the formation of the cluster takes place. 



the same time. The mean value of cr for each model is given 
in Col. 9 of Tabled Prestellar cores where the fit routine fails 
or where cr > 1 are not taken into account in our subsequent 
analysis. This concerns a wide variety of cores, there is no cor- 
relation with the final mass or the time of formation. However, 
they represent only a small fraction of the total number of ob- 
jects. The actual numbers of fitted cores are listed in Col. 8 of 
TablelU 

When interpreting the fit parameter logMo, the applied 
shift has to be taken into account. Thus, log Mjj,^^ = log(M())-7 
gives the real value of the peak accretion. This parameter is 
plotted for all protostellar cores and all models versus the re- 
spective final mass (Fig.|3j- A correlation with Mend is obvious. 
We apply a linear fit in the log-log diagram, which is indicated 
by the straight line. The fitted peak accretion rates show the 
same behaviour as the mean accretion rates (M). 

The parameter t indicates the time of the maximum of the 
accretion curve. It is plotted for all protostellar cores in Fig.|4] 
In almost all models t shows a coiTelation with the final mass. 
The parameter indicates how fast the gas falls in onto the core, 
therefore we expect it to be related to the local free-fall time 
and, thus, to the local density at the onset of collapse. It lies 
in the range 10"* < r < 10'' yr, which is less than the global 
free-fall time fjf . If we take an average value (r) ^ fff /3, this 
suggests an initial overdensity of almost a factor of ten in the 
collapsing regions. 



3.4. Class lifetimes and the effect of the turbulent 
medium 

We calculate the transition times from Class to Class I, as- 
sumed as described above. This gives the duration of Class 
phase for each protostar, the values are plotted versus the coiTe- 
sponding final masses in Fig.|5] The duration of Class phase 
increases with increasing final mass. Thus, a massive star is 
marked not only by a higher peak accretion rate but also by a 
longer time spent in Class phase. 

The mean accretion rates (M) of all individual protostel- 
lar cores of one model are averaged in four mass bins: < 
Mend/M© < 0.7 (binl), 0.7 < Mend/M© < 1.5 (bin2), 1.5 < 
Mend/M© < 3 (bin3), and Mend/M© > 3 (bin4). The values are 
given in Cols. 10 to 13 of Tabled Figure |6l shows the relation 
of the averaged mean mass accretion rate {M}mam to the Mach 
number for all models, split into mass bins and wave numbers, 
respectively. Three conclusions can be drawn from the figure: 
Firstly, there is a trend that (M)mean decreases with increasing 
IVIach number. That means that the mean accretion rate is lower, 
when the rms velocity dispersion (i.e. the turbulent Mach num- 
ber) of the medium is increased. The stronger support of the 
turbulent medium against gravitational collapse typically re- 
sults in a lower mass accretion rate. Secondly, (M)niean is higher 
for higher mass bins. This is consistent with the findings for 
the individual mean and maximum accretion rates ((M) and 
^max) discussed above. Finally, though, there is no coiTelation 
of (M)inean with the wavcnumbcr. Apparently the scale of the 
driving energy has no influence on the accretion rate. 
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4. Comparison with other theoretical models 

In the standard theory of isolated star formation (Shu ll977t . 
which takes a singular isothermal sphere as initial condition, 
the mass accretion rate is constant in time: M - 0.975 c^/G. 
Note that also the Larson-Penston solution (Larson 1969 
Penston '1969't, considering constant-density initial perturba- 
tions, gives a time-independent accretion rate, however, at 
a higher level of 47 c^/G. First numerical collapse calcula- 
tions have been reported by Bodenheimer & Sweigart fl968 1, 
Larson (1969 1, and Hunter J1977> . Models with more real- 
istic initial density profiles generally predict accretion rates 
that decline with time: The models of Hunter CI 977' Mmax 
» 36 Cj/G from numerical integrations of isothermal col- 
lapses), Foster & ChevaUer iTW^ M^rx ~ 47 c^/G from nu- 
merical hydrodynamic simulations), Tomisaka (1996' Mnax = 
(4-40) c^/G from numerical MHD models), Basu ( 1997 . Mm^x 
= 13 Cj/G from a semi-analytical model), Ogino et al. J 19991 
Mmax ~ (30 - 230) Cj/G from numerical hydrodynamic simu- 
lations), Masunaga & Inutsuka (2000, radiation hydrodynamic 
numerical codes), Whitworth & Ward-Thompson (120011 . and 
Motoyama & Yoshida (2003' M^ax ^ 42 /G) predict mass ac- 
cretion rates that peak shortly after the formation of the proto- 
star and decrease with time. This shows better agreement with 
observational data than constant accretion rates (see Sect. |5}- 
Our results display the same behaviour and our values of M^ax 
~ (3 - 50) Cj/G coincide quite well with those findings. In 
contrast, some models yield mass accretion rates that increase 
with time (McLaughHn & Pudritz 1997. Bonnell et al. .2001a 
Behrend & Maeder|200D- 

Theoretical models of star formation usually are scale-free. 
Thus, the results strongly depend on the adopted physical scal- 



ing, e.g. the choice of the initial density or temperature. The 
comparison of numerical values of accretion rates therefore re- 
quires some care. Again, in most cases the maximum accretion 
rates scale approximately as a few times Mj/rff . 

Whitworth & Ward-Thompson J2001> presented an analyt- 
ical model for protostellar collapse using a Plummer-like den- 
sity profile as initial condition. They successfully modelled the 
prestellar core L1544 in good agreement with observations. 
Their Mm!ix~ 8.1 x 10"^ M0yr"' coiTesponds quite well to 
our Gaussian collapse cases for the same stellar mass (Gl, 
G2). However, the accretion history of the collapsing Plummer 
sphere cannot be matched with our fit formula (|5}- The in- 
crease is steeper, while the decrease is slower compared to any 
of our models. The slow decline might be due to the fact that 
Whitworth & Ward-Thompson ( l20()Tt use a not-truncated, infi- 
nite density profile, while our models have finite sizes. A simi- 
lar model was used by Motoyama & Yoshida ( 2003|l who exam- 
ined the hypothesis that very high mass accretion rates exceed- 
ing 10""* M0yr"' require external triggering, as infeiTed from 
some observations. They find that the maximum accretion rate 
is proportional to the momentum given to the cloud core in their 
perturbed collapse model. A momentum of > 0.1 Mokms"' 
causes an accretion rate of > 10""* Mq yr"' . 

Smith (11999! I2000> presented a formula for the mass ac- 
cretion rate with a sharp exponential rise and a power law de- 
crease in time. This model provides an early peak in which 
M X! 10""* Mgyr"' for 10^ years and eventually becoming 
M X 10"' M0 yr ' for lO*" years. However, his formula (Eq. 6 
in Smith '2000 1 applies to our models only when choosing the 
parameters completely different than suggested, otherwise his 
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accretion curve has a more rapid increase but a slower decline 
than our models. 

Bonnell et al. (2001a) analysed competitive accretion in 
embedded stellar clusters by means of SPH simulations. They 
find that accretion in a cluster is highly non-uniform and that 
the accretion rate is higher for stars near the cluster centre. We 
do not see this in our results, likely because the protostars in 
our model are at different stages of evolution at a certain time, 
so this effect, if existent, is superposed by the strongly time- 
dependent variation in the accretion rate. Also the evolution of 
M with time in their models differs from our results: The mean 
accretion rate reported by Bonnell et al. (2001a), determined 
from all protostars in the cluster, increases with time until near 
the end of the simulation when the gas is significantly depleted. 
The difference is probably caused by the different assumptions 
(e.g. lack of turbulence, clustered potential). Indeed, the recent 
models of Bonnell et al. ( I2003t produce nearly constant accre- 
tion rates onto the most massive stars in the cluster The differ- 
ence might be due to the fact that the accretion rates are deter- 
mined by the accretion onto the cluster from outside, while in 
the models of Bonnell et al. (1200 lal all the mass was already 
in the cluster. 

Reid et al. ( I2002> used a logatropic equation of state as ba- 
sis for their hydrodynamical simulations of isolated star for- 
mation. Their accretion rate depends on the size of the core. 
It increases cubically and reaches the maximum when the 
expansion wave leaves the core, then it falls steadily. With 
the adopted scaling, M^ax is one to two orders of magnitude 
smaller than in our models. Consequently the whole accretion 
process lasts much longer, several 10^ years, which is in contra- 
diction to estimates of rapid star formation (Elmegreen 2000 
Hartmann et al. 120011 Hartmann l2003l Mac Low & Klessen 

Wuchterl & Tscharnuter ( I2003> find from models based 
on radiation hydrodynamics time varying accretion rates of a 
few 10"^ M0 yr"' for the phase < 0.8 Mend, which is about an 
order of magnitude lower than our values, especially for the 
Gaussian collapse or large-scale turbulence. The reason may 
be that protostellar cores in our models form by external com- 
pression before gravity takes over This results in enhanced ac- 
cretion rates relative to cores that begin contraction close to 
hydrostatic equilibrium. 

Hennebelle et al. ( I2003> performed numerical simulations 
where the collapse is triggered by a steady increase in the exter- 
nal pressure. M^ax is reached immediately after the formation 
of the central protostar (i.e. during Class phase), followed by 
a steady decrease to the Class 1 phase. The more rapid and the 
more prolonged the increase in external pressure, the higher is 
Mmax, ranging from 6.5 x 10"'' M© yr"' to 2.6 x 10"^ Mq yr"\ 
corresponding to ~ (4 - 16) /G. The qualitative behaviour of 
the accretion process does not differ much from our models, 
but the peak accretion rates are slightly smaller except in their 
models with the most rapid compression. 

5. Comparison with observations 

It is very difficult to measure mass accretion rates directly from 
observations (e.g. from inverse P Cygni profiles), instead they 



often have to be inferred indirectly based on the spectral en- 
ergy distributions (SEDs) of protostars or using outflow char- 
acteristics (Hartigan et al. 119951 Bontemps et al. I1996t . The 
correlation between accretion rates and outflow strength, how- 
ever, is still subject of strong debate (Wolf-Chase et al. l2003> . 
For Class objects typical mass accretion rates are estimated 
in the range 10"^ < MjUQyx-^ < 10"^ (Hartmann .19981 
Narayanan et al. [T998l Andre et al. [T999l CeccarelU et al. l2000l 
Jayawardhana et al.'200r, Di Francesco et al. 2001 ; Maret et al. 
|2^ Beuther et al. i2002aii2002bJ . The growth rate of Class I 
objects is believed to be about an order of magnitude smaller 
(Henriksen et al. 119971 Andre et al. 1200 1, with observational 
values between ~ 10"^ and ~ 5 x lO"*" M0yr"' (Brown & 
Chandler 1999] Greene & Lada l200l Boogert et al. l2002l 
Yokogawa et al. l2003l Young et al. l2003l . 

Bontemps et al. (I1996> studied the outflow activities in a 
sample of 45 low-mass embedded young stellar objects. They 
estimate that the observed decline of CO outflow momentum 
fluxes with time results from a decrease of the mass accre- 
tion rate from ~ 10"^ MQyr"' for the youngest Class pro- 
tostars to ~ 10"^ Mq yr"' for the most evolved Class I objects. 
Furthermore, they propose a simple exponential dependency 
of the accretion rate with time: M - (Mg,„/T)e"'^^ with ini- 
tial mass of the dense clump M*'jj^, and a characteristic time 
T w 9 X 10"^ yr. This is comparable to our Eq. (jSj- A similar 
exponential equation is also used by Myers et al. ( 1998 1, while 
Henriksen et al. ( 1997 i describe the accretion rate by an equa- 
tion that asymptotically approaches a power-law dependence at 
late times. 

Brown & Chandler ('1999V who determined an upper limit 
of M < (2 - 4) X 10"^ Mq yr"' for two Class I protostars in 
Taurus, also conclude that the accretion rate is not constant in 
time and likely much higher in the early phase. On the other 
hand, Hirano et al. (■2003J observed a dozen of deeply embed- 
ded young stellar objects of both Class and I type and derived 
the same mass accretion rates of (1 - 5) x 10"'' Mq yr"' for all 
of them. Unlike other authors, they argue that there is no sig- 
nificant difference in M between Class and Class I sources. 

The values given above correspond to the accretion rates 
derived for the model of gravoturbulent star formation dis- 
cussed here. They also decrease from 10"^ to lO"'* M0yr"' 
during Class phase to less than 10"^ Mq yr"' in later stages. 
However, the supposed transition between Class and Class I 
takes place still during the peak accretion phase. The accretion 
rates in our models typically do not decline significantly until 
about 80% of the final mass have been accreted (Fig.[0. This 
is unlike e.g. the model of Reid et al. (,2002). where M begins 
to fall off when about half of the mass of the core has been ac- 
creted. Given the uncertainties of the mass estimate for Class 
O/I transition we do not consider this a large discrepancy. 

According to observations. Class objects have an esti- 
mated Hfetime of ~ (1 - 3) x 10"^ yr (Andre et al. 2000 1. In our 
models this parameter is widely spread, ranging from ~ 10^ 
to > 10^ yr, but for a 1 M0 star it lies roughly in the range 
deduced from observations (see Fig.|5}- 
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6. Summary 

We have studied protostellar mass accretion rates from numeri- 
cal models of star formation based on gravoturbulent fragmen- 
tation. Twenty-four models covering a wide range of environ- 
mental conditions from low to high turbulent velocities and dif- 
ferent driving scales with a total number of 1325 protostellar 
cores have been investigated. Our main results may be sum- 
marised as follows: 

1. An order-of-magnitude estimate for mass accretion rates 
resulting from gravoturbulent fragmentation is given by 
M ^ Mj/Tff with Mj being the mean thermal Jeans mass 
and Tff the corresponding free-fall time. 

2. However, protostellar mass accretion is a highly time- 
variant process. It can be approximated by the empirical 
function logM(f) - logMo (e/r) f e"'^'". The peak accretion 
rate is reached during Class stage, shortly after the for- 
mation of the core; its value ranges between about 5 x 10"^ 
and 10""^ Mq yr"' . The maximum accretion rate is approxi- 
mately one order of magnitude higher than the constant rate 
predicted by the collapse of a classical singular isothermal 
sphere. 

3. Around the peak accretion phase the mass accretion rates 
are roughly constant. The mean accretion rates increase 
with increasing final mass. More massive stars have higher 
mass accretion rates and tend to form first. 

4. The same applies to the fitted peak accretion rates, which 
are also proportional to the final stellar mass. 

5. There is a similar correlation between the duration of 
Class phase (assuming that half of the final mass is ac- 
creted in this phase) and the final mass. 

6. (M)mean decrcascs with increasing Mach number of the tur- 
bulent environment, but is not correlated with the driving 
wavenumber. 

Our results agree well with many other models concerning the 
time evolution of the mass accretion process and the value of 
the peak accretion rate. In particular, the accretion rates from 
our models show an exponential decline, as it is also proposed 
by Bontemps et al. (.1996) . Myers et al. (il998j and Smith 
<1999lE50 0t. They also match observational findings like the 
supposed decline of the mass accretion rate from Class to 
Class 1 phase. We conclude that a theory of star formation based 
on gravoturbulent fragmentation of molecular clouds is an ad- 
equate approach to describe stellar birth in the Milky Way. 
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